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Biological systems fundamentally exist out of equilibrium in order to preserve organized struc- 
tures and processes. Many changing cellular conditions can be represented as transitions between 
nonequilibrium steady states, and organisms have an interest in minimizing the energy expended 
in such transitions. Using the Hatano-Sasa generalized dissipation, we extend a recently developed 
geometrical framework for determining optimal protocols so that it can be applied to systems driven 
from nonequilibrium steady states. We calculate optimal protocols for a colloidal particle dragged 
through solution by a translating optical trap with two controllable parameters, numerically ver- 
ify our optimal protocols, and offer specific experimental predictions. These optimal protocols may 
point to design principles for biological energy transduction systems and guide the design of artificial 
molecular machines. 



PACS numbers: 05.70.Ln, 02.40.-k,05.40.-a 

Introduction. Living systems are distinguished by 
their self-organization. Given the entropic driving force 
embodied in the second law of thermodynamics, creating 
and maintaining such organization requires staying far 
from equilibrium pQ, typically by coupling to nonequi- 
librium gradients. For example, ATP-driven molecular 
motors {e.g., kinesin) are forced away from equilibrium 
by cellular maintenance of a chemical potential differ- 
ence between ATP and ADP [2], and the rotary Fo-Fl 
ATP synthase operates out of equilbrium due to cellular 
maintenance of an electrochemical gradient across the in- 
ner mitochondrial membrane [3]. For constant ATP and 
ADP concentrations, or constant membrane potential, 
the dynamics of an ensemble of such molecular motors 
will approximate a nonequilibrium steady state (NESS). 
Thus, biological systems are often better characterized as 
nonequilibrium steady states rather than equilibrium sys- 
tems. Such NESS may change in response to changing en- 
vironmental conditions. Given that selective advantage 
may be incurred by energetically-efficient operation, evo- 
lution may have sculpted biological components to inter- 
act so as to reduce the energy wasted during transitions 
between NESS. Accordingly, optimizing such transitions 
can offer insights into the design principles of biological 
systems and guide the creation of synthetic molecular- 
scale machines. 

Though recent years have seen significant develop- 
ments in nonequilibrium statistical mechanics, a unify- 
ing set of principles governing nonequilibrium physics 
remains elusive. Much recent attention has focused on 



predicting optimal protocols to drive systems between 
equilibrium states with minimal expended work [41410J . 
In particular, Ref. [4] proposes a linear response frame- 
work for protocols that minimize the dissipation dur- 
ing nonequilibrium perturbations of microscopic systems. 
This idea is developed further in [5] where the utility 
of Riemannian geometry suggested in [4] is exploited to 
find explicit optimal protocols for a paradigmatic col- 
loidal particle model. 

In this article, we demonstrate that this geometric 
framework for obtaining optimal protocols may be ex- 
tended to NESS by optimizing the Hatano-Sasa dissipa- 
tion, sometimes referred to as the 'Y- value.' This mea- 
sure of dissipation plays a prominent role in a signifi- 
cant NESS fluctuation theorem that has been experimen- 
tally verified [T2]. The study of NESS promises greater 
mathematical tractability than arbitrary nonequilibrium 
phenomena, while significantly relaxing the restrictive as- 
sumptions of equilibrium physics. We apply this new op- 
timization framework to the paradigmatic model tested 
in [12] and analyzed in [13] with an eye towards ex- 
perimental tests. We calculate closed-form expressions 
for both the geodesic optimal protocol and the optimal 
straight-line protocol and test these protocols numeri- 
cally via a system of equations derived from the Fokker- 
Planck equation. Finally, we propose a regime of validity 
of our approximation based on this numerical work. By 
measuring the average work required to drive this system 
along either optimal or naive paths through control pa- 
rameter space, our results can be tested experimentally in 
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a straightforward way using existing experimental tech- 
niques. 

Derivation of the inverse diffusion tensor. We con- 
sider systems that relax to a stationary state distribu- 
tion p ss (x;X) when the set of control parameters A is 
held fixed. Here, x labels the microstates of the sys- 
tem. Stationary states include both equilibrium states 
(in which p ss is the familiar Boltzmann-Gibbs distribu- 
tion) and NESS. It will become apparent that, when the 
former case is considered, the inverse diffusion tensor de- 
rived in [5] coincides with the one presented below. 

The Hatano-Sasa dissipation arises in NESS transi- 
tions when the control parameters A are changed rapidly 
compared to the system's relaxation timescale. In some 
simple cases this corresponds to the system 'lagging' be- 
hind the changing control parameters. For transitions 
between equilibrium states this measure reduces to the 
standard dissipation governed by the Clausius inequal- 
ity |14| . In this paper we derive a framework, exact in 
the linear response regime |15j . for optimizing this dis- 
sipation during finite-time transitions between different 
nonequilibrium states. 

We define the ensemble average of the dissipation as 
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where 4>{x; A) = — In p ss (x; A) and r is the protocol dura- 
tion. During the driving process, the system's probability 
distribution over microstates fundamentally depends on 
the history of the control parameters A, which we denote 
by the control parameter protocol A. We assume the pro- 
tocol to be sufficiently smooth to be twicc-differentiable. 

Applying linear response theory [4j [5] [T5] and assuming 
that the protocol varies sufficiently slowly [4], we arrive 
at an expression for the dissipation 
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in terms of the control parameter velocities dX/dt and 
the inverse diffusion matrix C(A) with entries 
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The angle brackets (■ • • )a represent an average over noise 
followed by a stationary state average over initial con- 
ditions using the distribution p ss (x;X). Note that, if 
p as {x\ A) is the equilibrium distribution, the inverse dif- 
fusion tensor of [5] is recovered. 

In general, the detailed balance condition is violated 
in NESS and so the matrix (d\i(f)(t')d\i<l>{G)j, may be 
asymmetric. Eq. ^ shows the use of Eq. Q (specifi- 
cally its symmetric part) as a metric tensor is not pre- 
cluded. However, (d\i(j)(t')d\i(l>(Q)) x is not a covariance 
matrix and so a general proof of positive-definiteness is 
lacking [4 . These considerations do not affect the model 
considered here but future work will address this issue. 



The model system and its inverse diffusion tensor. 
We consider a particle with spatial coordinate x diffusing 
under Langevin dynamics subject to a one-dimensional 
harmonic potential, with equation of motion 
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x + n(t) - v(t) , 



for Gaussian white noise rj (t) satisfying 
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Here 7 is the Cartesian friction coefficient and x is the co- 
ordinate of the colloidal particle in the frame co-moving 
with the trap. The particle is initially in NESS due to 
constant trap velocity v. We take as our control parame- 
ters: the velocity v of the harmonic potential minimum, 
and the stiffness k of the trap (see Fig. 1 in Ref. |12|). 
The steady-state distribution is given by p ss (x;X) = 

f£ exp{-^ (kx + -fv) 2 } QHCd). The parameter space 
derivative of <fi is given by 
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In order to calculate the time correlation functions in 
Eq. pi), we solve Eq. Q for constant k and v, giving 
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Recalling that rj(t) is Gaussian noise, Eq. Q implies 
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Integrating over time yields the inverse diffusion tensor: 
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Optimal protocols. Though one can write down the 
geodesic equations for the metric [Eq. Q] in the (k, v) 
coordinate system, more insight is gained by finding a 
suitable change of coordinates. A direct calculation of 
this metric's Ricci scalar yields R = 0, demonstrating 
that the underlying geometry is Euclidean |16) . 



3 



The line element corresponding to the metric in Eq. ^ 
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To find the explicit isometry of this geometry into Eu- 
clidean space, we write the line element as 
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This suggests the coordinate transformation £ 



k ' 



so that 



ds 2 =f3 1 3 (de+d X 2 ) 
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In this coordinate system, geodesies are straight lines 
of constant speed. To find optimal protocols in (k, v) 
space, one simply transforms the coordinates of the 
endpoints into (£, x) space, connects these points by a 
straight line, and uses the inverse transformation to map 
the line onto a curve in (fc, v) space. This follows from 
the invariance of the geodesic equation |16j . Explicitly, 
the optimal protocol joining (fc^i^) and (kf,Vf) is 
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where T = — . Sample optimal protocols are pictured in 
Fig.|2] 

Optimal straight-line protocols. In the absence of any 
particular information about the system's dynamical 
properties, a naive control strategy would change the con- 
trol parameters at a constant rate, producing a straight 
line in control parameter space. The inverse diffusion ten- 
sor approximation [Eq. ^} provides a recipe for choosing 
both a potentially nonlinear path through control pa- 
rameter space, as well as a time-course along that path. 
The inverse diffusion tensor formalism can alternatively 
be used to optimize the time-course along a straight- 
line control parameter path. Such a protocol provides a 



benchmark against which an optimal protocol [Eq. ( 13 1] 



can be compared. For the model considered here, we will 
find that an optimal straight-line protocol can be sub- 
stantially better than the most naive (constant-speed) 
straight-line protocol. Furthermore, straight-line proto- 
cols are relatively straightforward to test experimentally. 

When k(t) is held fixed, a straightforward applica- 
tion of variational calculus demonstrates that a straight- 
line protocol in v(t) is exactly optimal and agrees with 
the predictions of the linear response approximation 



[Eq. pi]. In Ref. [T2], the average dissipation was mea- 
sured for three distinct experimental trials involving pro- 
tocols with constant k. As summarized in Fig. [T] the 
optimal protocol, namely the naive straight line in the 
case of constant k, shows significantly reduced dissipation 
compared with the protocols used in each experimental 
trial. However, in terms of testing the performance of the 
optimal protocols [Eq. ^M], kf ^ ki is the more general 



As in the case of finding globally optimal protocols, the 
problem of finding optimal straight line protocols simpli- 
fies dramatically in (£,x) coordinates. Using Eq. (12 1, 
we find 
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FIG. 1. Optimal protocols outperform constant-fc protocols 
tested in Ref. [T3]. Experiment 1 (left, red) used a quarter- 
sine wave protocol to vary the trap speed; Experiments 2 
(middle, blue) and 3 (right, orange) used an inverted three- 
quarters sine wave. Specifically, v(t) = 8.12 + 4.03 sin {wt/2r) 
for Experiment 1, v(t) — 9.93 — 3.63 sin (37rf/2r) for Experi- 
ment 2, and v(t) = 7.53 — 2.67 sin (3nt/2r) for Experiment 3. 
Here, velocity is measured in /im/s and r is the protocol du- 
ration measured in s. The average dissipation for these proto- 
cols (light color bar) and for the optimal protocols (solid color 
bar) were obtained numerically assuming fi- 1 = 4.6 pN nm 
(red), /3" 1 = 4.45 pN nm (blue), fi' 1 = 4.35 pN nm (or- 
ange) respectively. These effective temperatures were to 
chosen to give the best match between experiment and nu- 
merical calculation, and may differ from room temperature 
(/3 _1 = 4.14 pN nm) because of local heating by the optical 
trap [T7]. We predict a significant reduction in generalized 
dissipation for optimal protocol driving under the conditions 
of the three experiments described in Ref. [T5] . 
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The Euler-Lagrange equation implies 



which determines an implicit expression for x(0 : 
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The relation \ — 7V ^?fc determines an implicit expression 
for k(t), and hence for v(t). 

Computing dissipation numerically. We validate the 



optimality of the geodesies [Eq. ( 13 )] and compare with 



optimal straight-line protocols by calculating the average 
dissipation directly by integrating in time the Fokker- 
Planck equation [IB] . 
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In full generality, the mean dissipation as a functional of 
the protocol A(i) = (k(t),v(t)) is 



dt 



/? m 2 - 







2k 2 \ k 



■ Piv(x)a 



(19) 



Here angled brackets denote averages over the nonequi- 
librium probability density p{x 1 t). By integrating 
Eq. ( 18 ) against x and x 2 , we find a system of equations 
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We solve these equations numerically and compare the 
performance of optimal straight lines against geodesies 



[Eq. ( 13 )] and naive (constant-speed) straight-line proto- 
cols in Fig. [2] We selected endpoints and physical con- 
stants based on those used in the experiments of Ref. [H] ; 
in this regime the inverse diffusion tensor approximation 
produces small relative error. Though there is only a 
marginal difference in performance between the optimal 
straight-line protocol and the geodesic for both sets of 
endpoints, there is a substantial benefit in using either 
over the naive straight line protocol. 
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The inverse diffusion tensor arises naturally from the 
Fokker-Planck equation. If we neglect terms involving 
second- and higher-order temporal derivatives, we obtain 
an approximate solution to the Fokker-Planck system: 
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Substituting this into the expression for the mean dis- 
sipation [Eq. (19)], we recover Eq. (10 1. The argument 



above suggests that the emergence of the inverse diffu- 
sion tensor from the Fokker-Planck equation may follow 
from a perturbation expansion in small parameters [5]. 

Discussion. We have extended the geometric frame- 
work for calculating optimal protocols presented in [H [S] 
to systems relaxing to NESS. As energy-transducing bi- 
ological systems are more faithfully described by NESS 
than by equilibrium statistics, this brings recent theoret- 
ical developments significantly closer to the behavior of 
in vivo biological systems. 

Using a linear response approximation, we found the 
optimal dissipation for a model system of a colloidal par- 
ticle (initially in NESS) dragged through solution and 
subject to a time-dependent harmonic potential. We 
took as our control parameters the velocity and spring 
constant of the harmonic potential. As in [5 , tools from 
Riemannian geometry revealed a useful coordinate trans- 
formation which greatly simplified the construction of op- 
timal straight-line protocols as well as geodesic protocols. 
These optimal protocols were tested numerically and the 
small relative error in the dissipation approximation for 
experimentally relevant choices of parameters is encour- 
aging. 

Using the inverse diffusion tensor approximation in 
general allows us access to the full power of Riemannian 
geometry in calculating optimal protocols. However, ex- 
perimental tests can assess the range of validity of the ap- 
proximation. Our alternate derivation of the inverse dif- 
fusion tensor via a "derivative-truncation" expansion [5 
suggests a greater robustness of the approximation. 
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FIG. 2. Geodesies describe protocols that outperform 
naive (constant-speed) straight-line paths in parameter space. 
Geodesies between fixed pairs of points in the (k, u)-plane 
and accompanying straight-line protocols are pictured in (a). 
The filled circles represent points separated by equal times. 
The open circles correspond to the optimal parametrization 
along the respective straight path. All dissipations were cal- 
culated using the Fokker-Planck system, Eq. (20 1. Here, 
7 = 0.1 pN s//an and /3" 1 = 4.6 pN nm to approximate the 
experiments of Ref. |12J . The protocol duration is chosen to be 
t — 100 s to ensure that the relative error j 1 — {Y}^ pro * / (Y) a\ 
is less than 1.4% for all protocols. Protocol endpoints were 
selected for experimental accessibility |18] . The relative per- 
formance of naive straight-line, optimal straight-line, and 
geodesic protocols are summarized in (b). 



In this paper we provide concrete theoretical pre- 
dictions for experiments — specifically, we find that 
geodesies, optimal straight-lines, and naive straight-line 
protocols all are substantially more efficient than the pro- 
tocols tested in Trepagnier, et al. Moreover, we demon- 
strate that for simultaneous adjustment of k and v, opti- 
mal straight-line protocols can perform substantially bet- 
ter than naive (constant-speed) straight-line protocols. 
The necessary methodology and experimental apparatus 
are well-established [12] to not only verify our predictions 
but to push beyond the near-steady-state regime. 

Given the greater generality embodied by the exten- 
sion to NESS, and the accuracy of this approximation 
for a standard model system, optimal driving protocols 
derived in this framework may give insight into evolu- 
tionary optimization of biomolecular energy transduction 
and should provide operating principles for the design of 
artificial molecular machines. 

P. R. Z. thanks Christopher Jarzynski for a useful dis- 



cussion of this work, and Vladislav Belyy and Ahmet 
Yildiz for information regarding optical traps. M. R. 
D. gratefully acknowledges support from the McKnight 
Foundation, the Hellman Family Faculty Fund, the Mc- 
Donnell Foundation, and the Mary Elizabeth Rcnnie En- 
dowment for Epilepsy Research. M. R. D. and P. R. 
Z. were partly supported by the National Science Foun- 
dation through Grant No. IIS-1219199. D. A. S. was 
supported by NIGMS Systems Biology Center grant P50 
GM081879. 



[1 

[2 
[3 

[4: 

[5 

[6 

[7. 

[8 

[9 

[10 

[11 
[12 



* pzulkowski@berkeley.edu 
t l dav id.sivak@ucsf. edul 

* deweese@berkeley.edu 
E. Schrodinger, What is Life? (Cambridge University 
Press, Cambridge, 1992). 

J. Howard, Mechanics of Motor Proteins and the Cy- 
toskeleton (Sinauer, Sunderland, Massachusetts, 2001). 

B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, 
and P. Walter, Molecular Biology of the Cell (Garland 
Science, 2002). 

D. A. Sivak and G. E. Crooks, Phys. Rev. Lett. 108, 
190602 (2012). 

P. R. Zulkowski, D. A. Sivak, G. E. Crook s, and M. R. 
DeWeese, |Phys. Rev. E 86, 041148 (2012)| 
D. K. Shenfeld, H. Xu, M. P. Eastwood, R. O. Dror, and 
D. E. Shaw, Phys. Rev. E 80, 046705 (2009). 

D. C. Brody and D. W. Hook, J. Phys. A 42, 023001 
(2009). 

A. Gomez-Marin, T. Schmiedl, and U. Seifert, J. Chem. 
Phys. 129, 024114 (2008). 

T. Schmiedl and U. Seifert, Phys. Rev. Lett. 98, 108301 
(2007). 

E. Aurell, C. Mejfa-Monasterio, and P. Muratore- 
Ginanneschi, Phys. Rev. Lett. 106, 250601 (2011). 
T. Hatano and S. Sasa, Phys. Rev. Lett. 86, 3463 (2001). 
E. H. Trepagnier, C. Jarzynski, F. Ritort, G. E. Crooks, 

C. J. Bustamante, and J. Liphardt, |Proc. N atl. Acad. 



[13 

[14 
[15 
[16 
[17 



Sci. USA 101, 15038 (2004)] 
O. Mazonka and C. Jarzynski, "Exactly solvable model 
illustrating far-from-equilibrium predictions," (1999), 

arXiv:9912.121. 

M. Kardar, Statistical P hysics of Particles\ (Cambridge 
University Press, 2007). 

R. Zwanzig, Nonequilibrium statistical mechanics (Ox- 
ford University Press, New York, 2001). 
M. P. do Carmo, Riemannian geometry (Birkhauser, 
Boston, 1992). 

E. J. Peterman, F. Gittes, and C. F. Schmidt, Biophys- 



ical Journal 84, 1308 (2003) 



[18] K. Neuman and S. Block, Review of Scientific Instru- 
ments 75, 2787 (2004). 



